library(plyr)
library(estimatr)
library(emmeans)
library(tidyverse)
library(magrittr)
library(broom)
library(ggbeeswarm)

# load ipsos table

d0 <- "https://github.com/thomasjwood/full_fact/raw/main/ff_ipsos.rds" %>% 
  url %>% 
  gzcon %>% 
  readRDS

d1 <- "https://github.com/thomasjwood/full_fact/raw/main/ff_mt.rds" %>% 
  url %>% 
  gzcon  %>% 
  readRDS %>% 
  mutate(
    out_agree = out_agree %>% 
      mapvalues(
        c("Tend to disagree",
          "Tend to agree"),
        c("Disagree",
          "Agree")
      ) %>% 
      factor(
        c("Strongly disagree",
          "Disagree",
          "Neither agree nor disagree",
          "Agree",
          "Strongly agree")
      ),
    out_true = out_true %>% 
      factor(
        c("True",
          "Probably true",
          "Not sure",
          "Probably false",
          "False") %>% 
          rev
      ),
    agree_num = out_agree %>% 
      as.numeric,
    true_num = out_true %>% 
      as.numeric,
    cond = cond %>% 
      factor(
        c("cond_items",
          "cond_misinfo",
          "cond_corr")
      )
  )

d1$com_num <- d1$agree_num %>% 
  add(
    d1$true_num
  ) %>% 
  divide_by(2)

l1 <- list(
  c("vaccine_cancer",
    "pichetto_shanty",
    "safrica_latrines"),
  list(
    c("UKHQ7A_1",  "UKHQ7A_2"),
    c("ARGHQ7A_1", "ARGHQ7A_2"),
    c("ZAHQ3A_1",  "ZAHQ3A_2")
    )
  )

t1 <- l1 %>% 
  pmap_dfr(
    possibly(
      otherwise = NULL,
      function(i, j)
        
        d0 %>% 
        select(
          QCOUNTRY, Panellist_ID_Allocated, 
          any_of(j)
        ) %>% 
        mutate(
          issue = i
        ) %>% 
        set_colnames(
          c(
            "country", "caseid", "format_1", "format_2", "issue"
          )
        )  %>% 
        filter(
          (format_1 %>%
            is.na  |
          format_2 %>%
            is.na ) %>% 
            not
          )
    )
  ) %>% 
  modify_at(
    c(1, 3:4),
    as_factor
  )


t1$format_com <- case_when(
  
  # uk
  t1$issue %>% 
    equals(
      "vaccine_cancer"
    ) &
  t1$format_1 %>% 
    equals("Yes") ~ "Text",
  t1$issue %>% 
    equals(
      "vaccine_cancer"
    ) &
  t1$format_2 %>% 
    equals("Yes") ~ "Visual",
  
  # uk
  t1$issue %>% 
    equals(
      "pichetto_shanty"
    ) &
    t1$format_1 %>% 
    equals("Yes") ~ "Text",
  t1$issue %>% 
    equals(
      "pichetto_shanty"
    ) &
    t1$format_2 %>% 
    equals("Yes") ~ "Visual",
  
  # south africa
  t1$issue %>% 
    equals(
      "safrica_latrines"
    ) &
    t1$format_1 %>% 
    equals("Yes") ~ "Text",
  t1$issue %>% 
    equals(
      "safrica_latrines"
      ) &
    t1$format_2 %>% 
    equals("Yes") ~ "Video",
  
  TRUE ~ NA_character_
  )

d2 <- d1 %>% 
  full_join(
    t1 %>% 
      select(
        -matches("\\_\\d")
      )
  ) %>% 
  filter(
    issue %>%
      is_in(
        c("vaccine_cancer",
          "pichetto_shanty",
          "safrica_latrines")
      )
  ) %>%
  mutate(
    cond2 = str_c(
      cond,
      format_com %>% 
        is.na %>% 
        ifelse(
          "",
          str_c(
            "_",
            format_com
            )
          ) %>% 
        str_to_lower
      ) %>% 
      factor(
        
        c("cond_misinfo",
          "cond_items", 
          "cond_corr_text", "cond_corr_video", "cond_corr_visual")
      )
    )


t2 <- d2 %>% 
  filter(
    wave == "wave 1"
  ) %>% 
  group_by(
    country, wave, issue
  ) %>% 
  nest 



t3 <- 1:1000 %>% 
  map_dfr(
    function(b){
      
      t2$bd <- t2$data %>% 
        map(
          ~sample_frac(., 1, T)
        )
      
      t2$mods <- t2$bd %>% 
        map(
          ~lm_robust(
            com_num ~ cond2, 
            data = .x
          )
        )
      
      t2$emm <- t2$mods %>% 
        map2(
          t2$bd,
          ~emmeans(
            .x, 
            trt.vs.ctrl ~ cond2, 
            data = .y
          )
        )
      
      t2$edf <- t2$emm %>% 
        map(
          ~extract2(., 1) %>% 
            tidy
        )
      
      t2$cont <- t2$emm %>% 
        map(
          ~extract2(., 2) %>% 
            tidy
        )
      
      t2 %>% 
        ungroup %>% 
        select(
          country, issue, cont
        ) %>% 
        pmap_dfr(
          function(
            country, issue, cont
          )
            
            cont %>% 
            mutate(
              issue = issue,
              country = country
            )  %>% 
            slice(2:3)
        )
      }
    )
    
t3$iter <- 1:1000 %>% 
  rep(each = 6)

t4 <- 1:1000 %>% 
  map_dfr(
    function(b){
      
      t2$bd <- t2$data %>% 
        map(
          ~sample_frac(., 1, T)
        )
      
      t2$mods <- t2$bd %>% 
        map(
          ~lm_robust(
            com_num ~ cond2, 
            data = .x
          )
        )
      
      t2$emm <- t2$mods %>% 
        map2(
          t2$bd,
          ~emmeans(
            .x, 
            trt.vs.ctrl ~ cond2, 
            data = .y
          ) %>% 
            pairs
        )
      
      t2$cont <- t2$emm %>% 
        map(
          ~extract2(., "contrasts") %>% 
            as.data.frame %>% 
            slice(3)
        )
      
      t2 %>% 
        ungroup %>% 
        select(
          country, issue, cont
        ) %>% 
        pmap_dfr(
          function(
            country, issue, cont
          )
            
            cont %>% 
            mutate(
              issue = issue,
              country = country
            )
        )
    }
  ) 

t4$iter <- 1:1000


t5 <- t3 %>% 
  mutate(
    horiz = "correction_effects"
  ) %>% 
  bind_rows(
    t4 %>% 
      mutate(
        horiz = "Difference in correction effects",
        contrast = t3$contrast[[1]]
      ) %>% 
      rename(
        std.error = SE,
        adj.p.value = p.value,
        statistic = t.ratio
      ) 
  )

t5$issue %<>% 
  mapvalues(
    c( "pichetto_shanty",
       "safrica_latrines",
       "vaccine_cancer"),
    
    c("4k shanty towns in Buenos Aires",
      "100s children drown in pit latrines",
      "Vaccines might cause cancer")
  )


t5$horiz %<>% 
  mapvalues(
    c("correction_effects", 
      "Difference in correction effects"),
    c("Correction effects", 
      "Difference in correction effects")
    ) %>% 
  factor(
    c("Correction effects", 
      "Difference in correction effects")
    )

t5$contrast %<>% 
  mapvalues(
    c("cond_corr_text - cond_misinfo",
      "cond_corr_visual - cond_misinfo", 
      "cond_corr_video - cond_misinfo"),
    c("Correction (Text) - Misinformation",
      "Correction (Visual) - Misinformation", 
      "Correction (Video) - Misinformation")
    ) %>% 
  factor(
    c("Correction (Text) - Misinformation",
      "Correction (Visual) - Misinformation", 
      "Correction (Video) - Misinformation")
    )

t5$issue %<>% 
  fct_reorder(
    t5$estimate, .fun = function(i) mean(abs(i)), .desc = T
  )

t5l <- t5 %>% 
  group_by(
    contrast, issue, horiz
  ) %>% 
  summarize(
    mu = estimate %>% mean
  ) %>% 
  mutate(
    lab = mu %>% 
      round(2)
  )

t5l$lab %<>% 
  str_extract("\\d") %>% 
  equals("0") %>% 
  ifelse(
    t5l$lab %>% 
      str_replace("0.", "."),
    t5l$lab
  )

t5l$lab %<>% 
  mapvalues(
    "-.3",
    "-.30"
    )

t5 %>% 
  ggplot() +
  geom_vline(
    xintercept = 0,
    size = .25,
    linetype = "dashed"
  ) +
  geom_quasirandom(
    aes(estimate, contrast),
    alpha = .05,
    groupOnX = F
    ) +
  geom_point(
    aes(
      mu, contrast, 
    ),
    shape = 21,
    size = 6.5,
    fill = "white",
    data = t5l
  ) +
  geom_text(
    aes(
      mu, contrast, label = lab 
    ),
    size = 1.75,
    data = t5l
  ) +
  facet_grid(
    issue ~ horiz, 
    scales = "free", 
    space = "free_x", 
    labeller = label_wrap_gen(width = 20) 
  ) +
  scale_y_discrete(
    breaks = t5$contrast %>% 
      levels,
    labels = t5$contrast %>% 
      levels %>% 
      str_replace(
        " - ",
        "-\n"
      )
  ) +
  labs(
    x = "",
    y = "",
    color = ""
  ) +
  theme(
    panel.grid = element_blank(),
    strip.text.x = element_text(
      angle = 0
    ),
    strip.text.y = element_text(
      angle = 0
    ),
    legend.position = "none",
    plot.margin = margin(.2, .2, .2 , 1, "cm") 
  )


